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Abstract 

We present an example of a highly connected closed network of 
servers, where the time correlations do not vanish in the infinite vol- 
ume limit. The limiting interacting particle system behaves in a peri- 
odic manner. This phenomenon is similar to the continuous symmetry 
breaking at low temperatures in statistical mechanics, with the aver- 
age load playing the role of the inverse temperature. 
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1 Introduction 

1.1 Interacting particle systems with long range mem- 
ory 

The theory of phase transitions, among many results, substantiates the pos- 
sibility of constructing reliable systems from non-reliable elements. As an 
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example, consider the infinite volume stochastic Ising model at low tempera- 
ture T in dimension > 2, see [L]. It is well known, that if we start this system 
from the configuration of all pluses, then the evolution under Glauber dy- 
namics has the property that the fraction of plus spins at any time exceeds 
i+m (T) ^ w ] 1 i c j 1 j g D igg er than | for T < T cr . (Here m* (T) is the spontaneous 
magnetization.) On the other hand, if we consider finite volume Ising model 
(with empty boundary condition, say), then this property does not hold, and 
the system, started from the all plus state, will be found in the state with the 
majority of the spins to be minuses at some later (random) times. Therefore, 
the infinite system can remember, to some extent, its initial state, while the 
finite system can not. 

There are many other examples of this kind, that belong to the theory 
of interacting particle systems, such as voter model, contact model, etc. In 
all these examples we see systems, that are capable of "remembering" their 
initial state for arbitrary long times. 

In the present paper we are constructing a particle system that "remem- 
bers its initial phase" . The rough analogy can be described as follows. Imag- 
ine a brownian particle (p (t) , with a unit drift, which lives on a circle. Sup- 
pose the initial phase p (0) = 0. Then the mean phase <p it) = tmod (2tt) , but 
with time we know the phase p (t) less and less precisely, since its variance 
grows, and in the limit t — > oo the distribution of p (t) tends to the uniform 
one. However, one can combine infinitely many such particles by introduc- 
ing suitable interaction between them in such a way that the memory of the 
initial phase does not vanish and persists in time. Namely, one has to put 
such particles at the sites of Z 3 and to introduce the attractive interaction 
between them. If the initial state is chosen to be coherent, then the phase of 
every particle will grow linearly, while its variance will stay bounded. 

This is roughly what we do in the present paper. We consider a network 
of simple servers which are processing messages. Since the service time of 
every message is random, in the course of time each single server loses the 
memory of its initial state. So, in particular, the network of non-interacting 
servers, started in the same state, becomes de-synchronized after a finite time. 
However, if one introduces certain natural interconnection between servers, 
then it can happen that they are staying synchronized after an arbitrary long 
time, thus breaking some generally believed properties of large networks. We 
have to add here that such a phenomenon is possible only if the mean number 
of particles per server is high enough; otherwise the infinite network becomes 
de-synchronized, no matter which kind of interaction between servers takes 
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place. So the parameter of the mean number of particles per server, called 
hereafter the load, plays the same role as the temperature in the statistical 
mechanics. 

In other words, the transition we describe happens due to the fact that 
at low load the behavior of our system is governed by the fixed point of the 
underlying dynamical system, while at high load the dominant role is played 
by its periodic attractor. A similar phenomenon was described by Hepp and 
Lieb in [HE] . 

Below we present the simplest example of the above behavior. But we 
believe that the phenomenon we describe is fairly general. Its origin lies in 
the fact that any large network of the general type possesses some kind of the 
continuous symmetry in the infinite limit, and it is breaking of that symmetry 
at high load that causes the long-range order behavior of the network. In 
our case this is the rotation symmetry, corresponding to the periodic orbit 
of the limiting dynamical system. 

1.2 Information networks and their collective behavior 

Now we will describe one pattern of behavior of certain large networks, which 
was assumed to be universal. It is known under the name of Poisson Hy- 
pothesis. 

The Poisson Hypothesis is a device to predict the behavior of large queu- 
ing networks. It was formulated first by L. Kleinrock in [K] . and concerns 
the following situation. Suppose we have a large network of servers, through 
which many customers are traveling, being served at different nodes of the 
network. If the node is busy, the customers wait in the queue. Customers 
are entering into the network from the outside via some nodes, and these 
external flows of customers are Poissonian, with constant rates. The service 
time at each node is random, depending on the node, and the customer. The 
PH prediction about the (long-time, large-size) behavior of the network is 
the following: 

• consider the total flow T of customers to a given node M . Then T is 
approximately equal to a Poisson flow, V, with a time dependent rate 
function A^v (T) . 

• The exit flow from M - not Poissonian in general! - has a rate function 
7_/V (T) , which is smoother than Ajv (T) (due to averaging, taking place 
at the node Af) . 
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• As a result, the flows (T) at various nodes M should go to a constant 
limits Xj\f ~ 7f Jq X (t) dt, as T — ► oo, the flows to different nodes being 
almost independent. 

• The above convergence is uniform in the size of the network. 

Note that the distributions of the service times at the nodes of the network 
can be arbitrary, so PH deals with quite a general situation. The range of 
validity of PH is supposed to be the class of networks where the internal 
flow to every node AT is a union of flows from many other nodes, and each 
one of these flows constitutes only a small fraction of the total flow to N '. If 
true, PH provides one with means to make easy computations of quantities 
of importance in network design. 

The rationale behind this conjectured behavior is natural: since the inflow 
is a sum of many small inputs, it is approximately Poissonian. And due to 
the randomness of the service time the outflow from each node should be 
"smoother" than the total inflow to this node. (This statement was proven 
in [RShVj under quite general conditions.) In particular, the variation of the 
latter should be smaller than that of the former, and so all the flows should 
converge to corresponding constant values. 

In the paper [RSh] the Poisson Hypothesis is proven for simple networks 
in the infinite volume limit, under some natural conditions. For systems with 
constant service times it was proven earlier in [St lj . 

The purpose of the present paper is to construct a network that satisfies 
the above assumption - that the flow to every given node is an "infinite" sum 
of "infinitesimally small" flows from other nodes - but has coherent states. 
That means that the states of the servers are evolving in a synchronous 
manner, and the "phase" of a given server behaves (in the thermodynamic 
limit - i.e. in the limit of infinite network) as a periodic non-random function, 
the same for different servers. 

We have to stress that our network exhibits these coherent states only in 
the regime when the average number N of the customers per server - called 
in what follows the load - is large. For low load we expect the convergence to 
the unique stationary state. This "high temperature" kind of behavior will 
be the subject of the forthcoming work. 

Our network Vqo is constructed from infinitely many elementary "tri- 
angular" networks V (described below, in Section 12.11) . A single triangle 
network V = Vi with iV customers is just a Markov continuous time ergodic 
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jump process with finitely many states. As N becomes large, this Markov 
process tends (in the appropriate "Euler" limit) to a (5-dimensional) dynam- 
ical system A, possessing a periodic trajectory C, which turns out to be a 
stable (local) attractor. The coordinate <p parameterizing that attractor C 
is the "phase" alluded to in the previous subsection. The combined network 
corresponds in the same sense to the coupled family Aqo of dynamical sys- 
tems A. We establish the synchronization property of that coupled family 
A^ and that allows us to construct coherent states of the network Vcx>. 

The networks Vm, composed of M triangle networks V, are ergodic. 
Their evolution is given by irreducible finite state Markov processes with 
continuous time. Let n M be the invariant measure of the process Vm- As 
M — > oo, the sequence of Markov processes Vm converges weakly on finite 
time intervals to a certain limiting (non-linear Markov) process Voo- By the 
theorem of Khasminsky - see Theorem 1.2.14 in [L] any accumulation point 
of the sequence ttm is a stationary measure of Voo- The special measure Xoo, 
describing "the Poisson Hypothesis behavior", is also a stationary measure 
of Voo- If Xoo is a global attractor of Voo, then, of course, the Poisson 
Hypothesis holds. The proof of the Poisson Hypothesis in |RSh] was based 
on this argument. The existence of an accumulation point of the sequence ttm 
that differs from Xoo would be the strongest counterexample to the Poisson 
Hypothesis. This problem will be addressed in forthcoming papers. Here we 
prove a weaker statement that Xoo is not a global attractor for Voo- 

In [RStJ Rybko and Stolyar observed that the condition that the work- 
load at every node of a multiclass open queueing network is less than 1 is 
not sufficient for the network to be ergodic. In connection with this, they 
introduced a new approach to the analysis of ergodicity of networks, which 
reduces the problem to the question of stability of the associated fluid mod- 
els. It was shown by them that the two-node priority network, considered 
in jRSt], is ergodic if and only if for every initial state of the corresponding 
fluid model the total amount of fluid vanishes eventually. This approach was 
further developed by Dai [Dj, Stolyar |St2j . and Puhalsky and Rybko [PR] , 
who proved that stability of the fluid model is necessary and sufficient for 
ergodicity of a certain class of general networks. Interesting instances of non- 
ergodic queueing networks with mean load being smaller than the capacity, 
where considered by Bramson [Bit IB2j . Our construction will be based on 
the following open network introduced by Rybko and Stolyar (RS-network) 
in [KStj . 

This queuing network with four types of customers is represented by the 
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following 4-dimensional Markov process. Customers arrive to the network 
according to Poisson inflows of constant rate A. The network consists of two 
nodes - A and B. All the service times are exponential, hence the network 
is defined by the rates, the evolution of types of the customers and the 
priorities. The customer of type A (respectively, B) arrives to the node 
A (respectively, B). The customer A is served with the rate 7a, then is 
sent to B, with type AB. There he is served with the rate 'Jab and leaves 
the network. Symmetrically, jb = 7a, and jba = Jab- Each customer 
AB is served before all the customers B, while each customer BA is served 
before all the customers A. The nominal workload at nodes A and B equals 
P = ^(1a + Tra)- The serv i ce rates satisfy the conditions jab < 2A and 
p < 1. It is proved in [RStj that for certain values of the parameters the 
resulting Markov process is transient. The fluid limit (or the Euler limit) of 
this network evolves in the following non-trivial manner: each node is empty 
during a positive fraction of time, but at other moments it is non-empty, 
and, moreover, the total amount of the fluid in the network tends linearly to 
infinity. 

The rest of the paper is organized as follows. In Section 2 we define our 
networks VjJJ-. Here M is the size of the network and N is the load per node. 
We formulate the preliminary version of our Main Result. In Section 3 we 
study the limiting network, V^, and prove the convergence Vjy- — > V^. In 
Section 4 we introduce the fluid networks, A M , which are coupled dynamical 
systems, and their limit, Aoo, which turns out to be a non-linear dynamical 
system, in the sense made precise in this Section. In particular, we show 
that Aqo is not ergodic. In the next Section 5 we prove the convergence of 
the Non-Linear Markov Process to its Euler fluid limit, Aqo, as N —>■ 00. 
The last Section 6 contains the formulation and the proof of our main result, 
Theorem [T7J 

To save on notation, we consider throughout this paper the simplest el- 
ementary symmetric model, depending on 3 parameters. We stress the fact 
that this (discrete) symmetry is not essential in our case, and our results are 
valid for any small 6D-perturbation of our model. 
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2 Mean-field network and its limit 

2.1 Basic network 

We will consider the following 5-dimensional Markov process, V N . It de- 
scribes a closed queuing network with N customers. It consists of three 
nodes: O, A and B, through which the customers go. All the service times 
are exponential, so we only need to specify the rates, the evolution of types of 
the customers and the priorities. To simplify the presentation we will make 
a specific choice of these rates. The node O serves all the customers on the 
FIFO basis, with the rate 70 = 3. After being served, the customer goes to 
the node A or to B, choosing one of them with probability |. If he arrives 
to A, he gets the type A, otherwise B. The customer A is served with the 
rate 7^ = 10, then is sent to B, with type AB. There he is served with the 
rate jab = 2 and goes back to O. Symmetrically, 75 = 10, and ^ba = 2. 
Each customer AB is served before all the customers B, and each customer 
BA is served before all the customers A. More precisely, if an AB customer 
arrives to the B node, while the node is serving some B customer, his service 
is stopped and is resumed only at the moment when the service of all AB 
customers is over. 



O 
O 

o 




The elementary network. 



Of course, the above choice of the rates is not the only possible. Any 
other choice would be as good, provided the corresponding fluid network, 
which can be associated to our queuing network, has some specific property 
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- namely, we need this fluid network to have a cyclic regime. The fluid 
network will be described in details in the Section ET21 below. 

2.2 M coupled processes 

Let be the Markov process, obtained from M copies of V^, intercon- 
nected in the mean-field manner. We take the total number of customers to 
be NM. 

The mean field network is defined as follows. Each node Oi, % = 1, M, 
is connected to all of the nodes Aj, Bj, j = 1,...,M, and each customer, 
leaving the node Oi, goes to each of the 2M nodes Aj, Bj with the same 
probability The rate of leaving the node Oi is the same, as above, i.e. 
equals to 70 = 3. In a similar way, the A customers of every node A\ are 
exiting it with the rate jab = 10, and then choose one of the Bj nodes with 
probability j?, and so on. The priorities are kept the same: if the node A\, 
say, is in the state with customers of both kinds - A and BA - present, then 
the BA customers are served first, with no delay. 




Two coupled processes. 



The configuration of the process is given by the number of customers of 
each type at each of the 3M nodes, that is by an integer point in (1R 5M ) + . 
Due to the mean-field symmetry, we can factor the set of configurations by 
the product of permutation groups Sm x Sm x Smi an d still have the Markov 
process. The orbit of the symmetry group corresponds to a collection of M 3 
integer points %i G (M 5 ) + , some of which may coincide. 
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It is convenient for us to index these configurations by the atomic mea- 
sures, 

M 3 



i=l 

In fact, they belong to the set M ^(-^Z 5 ) + j . Note that every such measure 
\i factor into a product 

H = (no, HA, Hb) = Ho x HA x Hb = ^6 [h] x ^-a [h] x ^-b [h] 

of probability measures on M 1 = {xo} , resp. IR 2 = {xa,%ba} and M 2 = 
{xb,%ab} ■ Here we denote by n*-s the various projections (or marginals). 
We have /io = ^i=i^' f° r some ( n °t necessarily distinct) Xi 6 Z 1 , i = 

1,...,M, likewise ^ = Ms = if E 2 =i%> e Z 2 . We will 

N N 

denote the set of all such measures by Mm- The state v M of our Markov 
process is then an element from M (Mm) , i-e. a measure on the measure 
space. Among these there are configurations of the process Vm, namely, the 
5- measures S m , with m G Mm, so we can define in this way the embedding 
Mm C M (M m) , However, even if the initial state v M (0) of Vm happen 
to be such a measure 5 m , i.e. v M (0) G A4m, then at any positive t we have 
only that z/ M (t) G M (Mm) , while in general v M it) ^ A4m- 

For the future use we will write down the rates of the factor-process. Let 
v be some measure of the form ([1]) , while v' be the measure obtained from 
v after a single jump of the initial process. For example, let us consider the 
case when the jump in question is of AB — > O type, from S-type server to 
O server (with the rate Jab)- That means that for some unique well-defined 
(by the pair v,v') elements x§ = (xb,%ab) G (^Z 2 ) , i G G (^Z 1 ) we 
have: ^ ^ 

t/ = v (x B ) - jj, v' (x ) = v (x ) - jj. (2) 
Of course, for another pair: x^ = (xb, xab — 1) , %o — x o + 1> we have 

v' (X B ) = V (x B ) + Jj, V' (X ) = V (X ) + "^r, (3) 

while at all other locations the two measures are the same. Since there 
are Mv (x§) locations where the jump could originate, and the fraction of 
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sites with the desirable outcome is v (xo) , we have for the rate c (v, v') = 
cab iy,v') the expression 



c (v, v') = jabMv (xg) v (x ) ■ 



If the measures v, v' are not related by (J2j) — ([3]) , then the rate cab {v, v') = 0. 
We will keep the notation V m for the factor-process. 

2.3 M — > oo limit: Non-linear Markov Process 

Suppose that a sequence of initial states v M (0) G M. (A4m) of the Markov 
processes Vm is given, which satisfy vm (0) = S mM , with uim G .Mm, and 
moreover the weak limit v = \im.M-*oo f^M exists. Then the weak limits 
v (t) = liniAf^oo vm {t) exist for every t, and, moreover, for every t we have 
v (t) G M.. This is the Non-Linear Markov Process, NLMP, Vqo. The process 
is called Non-Linear since the transition mechanism to evolve from a given 
configuration depends not only on that configuration, but also on the measure 
from which this configuration was drawn. Such processes were introduced 
in [MB lM2], see also [RSh] . The above limiting NLMP-s depend on the 
parameter N, which is the number of clients per basic queuing network. We 
want to study the dependence on N, so we explicitly (re)introduce the index 
N in our notation. Thus, u N (t) refers to the states of the process V^. 

We will describe the limiting NLMP in the next Section [31 Now we can 
formulate the preliminary version of our main result. 

Theorem 1 Consider the Non-Linear Markov Process Voo, started from the 
measure Uq , which is close enough to the atomic measure with the single 
atom at vector X (A, N) G (-^Z 5 ) + , having coordinate xa — 1 and all other 
coordinates zero. "Close enough" here means that for some e > small 
enough we have Pkrov [Vq, o~jc(A,N)) < £ - Suppose additionally that the a- 
exponential moment of the measure Vq is less than a certain quantity E; 
a = a (e) , E = E (e) . Then the measure vf does not converge to any limit 
as t — > oo, provided N is large enough. 

More precisely, there exists a sequence of times t' k — > oo, such that 




., N)) is a neighborhood of X (A, N) 
oo. At the same time, there exists 
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another sequence t'l — > oo, for which v^, [Z7jv (X (B, N))] > 1 — 8^. In words, 
the measure Vt exhibits oscillations. 

Accordingly, the states of finite size networks, (Hkf)* > exhibit oscillations 
for long times, before going to their limits. The duration of the oscillation 
regime diverges with M. Different components of (v M ) t are oscillating almost 
coherently, for large M. 

3 The convergence —> : application of 
the Trotter-Kurtz theorem. 

Here we prove the convergence of the Markov processes V M to the Non-Linear 
Markov process V^. We will do that by writing down their generators Am 
and A, and by subsequent application of the Trotter-Kurtz theorem (Prop. 
1.3.3 in [EKJ), which we formulate now. 

Let Am, A : X — > X are (unbounded) operators on the Banach space X, 
and X C X is a dense subspace, belonging to the domains of definition of all 
Am-s and A. The following two conditions are sufficient for the convergence 
of the semigroups exp {tA M } — > exp {tA} on X as M — > oo : 

1. VV> G X we have A M (ip) — > A (ip) as M — > oo; 

2. there exists a dense subspace X\ C X , such that G X x we have 
exp {tA} (ip) G X . Such subspace X is called a core of A. 

3.1 Equation for the evolution V^. 

Here we study the limiting process V^. We write down its generator, and 
we exhibit its core. 

Let Vt be the evolution of the measure under V^. To find it we have to 
specify the initial measure Uq and then to solve the Cauchy problem for the 
differential equation, which equation we will write now. 

To do it we first introduce the (Poisson) rates 
A (t) = (Xo (t) , Xa {t) , A_b (t) , Xab (t) ,XBA(t)), corresponding to the state 
v t ■ 

X a (t) = 7a Ut ' fOT a = > AB > BA ' W 

x:x a >0 
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Aa {t) =1a ^ (x) , Ab (£) = 7b Y Vt ^ ■ ( 5 ) 

x:x A >0,x BA =0 x:x b >0,x ab =0 

We also introduce the 5D vectors A a to be the basis vectors of the lattice 
Z 5 . Then 



dt 



\a=O.A,B,AB,BA 



ft 0*0 ( 2^ la(l-Sx a )+lA(l-Sx A )Sx BA +lB(l-Sx B )Sx A1 
\a=0,AB,BA 

+ u t (x- A ) (1 - S wo ) (X AB (t) + Xba it)) (6) 
+ £ , t (x-A )(l-0^ 

a=A,B 

+ i/ t (x - A AB ) (1 - 8 XAB ) \ A (t) +v t (x- A BA ) (1 - 5 XBA ) X B (*) 
+ Y, vt{x + A 

a) 7a 

a=0,AB,BA 

+ u t (x + A A ) 1a&x ba + i> t (x + A B ) j B Sx AB ■ 

This is the value of the function A(p x , where A is the generator of the 
Markov semigroup S t = exp {tA} of the process V^, and the function <p x , 
which on every measure v G M. (Z 5+ ) takes value v (x) , computed at the 
point v t . As we will show below, the system (jSJ) has a unique solution. 

For reasons which will be explained later, it will be more convenient for us 
to use another basis in the space of functions on measures. Namely, for every 
v G M. (Z 5+ ) and every x G Z 5+ we define the function u (x) = J2 y > x v (v) > 
where the summation goes over all sites y such that all the coordinates of 
the difference y — x are non-negative. Then the functions A a (t) (see (jlj)) are 
given in the new variables as 

Aa (t) = laU t (A a ) , 

while the action of the generator A on the function u is given by the following 
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(simpler) equation: 

du t (x) 



a=0,AB,BA 

- (u t (x) -ut(x + A A ) -u t {x + A AB ) + u t {x + A A + A AB ))^ A (1 - 8 XA ) 5 XBA 

- (m (x) -u t {x + A B ) -u t {x + A BA ) + u t {x + A B + A BA )) ^ B (1 - 5 XB ) 5 XAB 
+ (u t (x - A ) - u t (x)) (1 - S Xo ) (j AB u t (A A b) + iBAUt (A B a)) 

(7) 

+ Y, Mx-A a )-u t (x))(l-8 Xa )^^ 

a=A,B 

+ (u t (x - A AB ) - u t (x)) (1 - 8 XAB ) j A u t (A A +A BA ) 
+ (u t (x - A BA ) - u t (x)) (1 - 5 XBA ) 'jsut (A B +A AB ) . 

The first three lines correspond to the second line of the equation ([6]) , while 
the last four - to the lines 3-5; the remaining lines of it disappear from 
the equations for u. The advantage of (J7j) over (jHJ) is that the equation for 
dut }^ contains only ut (y)-s with y-s in some finite set Y (x) , and moreover 
max, \Y(x)\= 20. 

Of course, the coordinates u (-)-s on Ai (Z 5+ ) are not independent. There 
are two kinds of relations between them: 

1. every value v{x) equals to L x (u) , where L is a certain linear form, 
depending on u (y) with y-s having form y = x + Y^ e a A a , e a = 0, 1; we 
need that for all x L x (u) > 0; 

2. 

lim u (x) = 0. (8) 

For technical reasons we will extend the action of our Markov semigroup 
to the space Ai (K) of measures on the compactification K of the lattice Z 5+ , 
where 

K= {Z+ + oo} 5 . 

The functions {u (x) ,x G Z 5+ } on Ai (IK) also play the role of coordinates 
there, provided that the relation (jSJ) is dropped. The evolution of the mea- 
sures is given by the same set of equations (jTj) . 

We supply Ai (K) with the topology of weak convergence. (We repeat for 
clarity that the subset Ai (Z 5+ ) C Ai (K) is invariant under our semigroup.) 
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Let C° = C {M. (K)) be the space of functions on M. (K) , continuous with 
respect to this topology. 

Theorem 2 The semigroup St acts on the space C = C (Ai (IK)) of continu- 
ous functions on M. (K) , and is strongly continuous and contracting. 

Proof. 1. Let us show the existence of the solutions to ([7]) . The equations 
((7j) are describing the evolution of the "closed" system. Consider now the 
corresponding "open" system, defined by the (arbitrary) rates A = A a (t) 
of the Poisson inflows and the initial state Uq. It evolves according to the 
equations 



< ^M = ~ (ut(x)-ut(x + A a )) la (l-5 Xa ) 

' a=0,AB,BA 

- (u t (x) -u t (x + A A ) -u t (x + A AB ) + u t {x + A A + A AB ))^ A (1 - 5 XA ) 5 XBA 

- (u t (x) -u t (x + A B ) -u t (x + A BA ) + u t [x + A B + A BA )) ^ B (1 - 8 XB ) 8 XAB 
+ (u t (x - A ) - u t (x)) (1 - 5 xo ) (X AB (t) + X BA (t)) 

+ (ut(x-A a )-u t (x))(l-5 Xa )^- 

a=A,B 

+ (u t (x - A AB ) - u t (x)) (1 - 5 XAB ) X A (t) 
+ (u t (x - A BA ) -u t (x))(l- 8 XBA ) \ b (t) . 

The corresponding exit rates b a are given by the natural relations 

bl it) = la u t (A a ) . 

Consider the function d x (t) : 

4 it) = b\ B it) + b\ A it) , 
d \ it) = 4 (t) = hi (t) , 

d\ B it) = b\ it) , d\ A (t) = b\ it) . 

The closed system is a fixed point of the map A d x , i.e. a solution of the 
equation 

A = J\ 
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To see the existence of a fixed point, let us introduce the functions A, B A and 

M0= f\(t)dt, B A (t) = fb\(t)(t)dt, D A (t) = fdl{t)dt, 
Jo Jo JO 

and the corresponding mapping A D A . The functions A., B A and D A are 
monotone continuous, moreover, the functions B A and D A have uniformly 
bounded derivatives. Let c be the upper bound for these derivatives, and 
CC be the space of all continuous monotone 5D vector-functions on [0, T] , 
vanishing at zero, with the derivatives bounded by c once they exist. Then 
<£ is compact and convex, therefore the map ^ uo : € — > € has at least one 
fixed point. 

2. We now will show that for every uq the map ty Uo is a contraction; that 
will imply the uniqueness of the solution. Without loss of generality we can 
assume that T is small. Informally, the contraction takes place because the 
exit rates (t) for t G [0, T] with T small depend mainly on the initial state 
u : the new clients, arriving during the time [0, T] have no chance to be 
served before T, if there were clients already waiting. Therefore the "worst" 
case for us is when the initial state u Q is the measure 5q, having a unit atom 
at G Z 5+ . 

So let Xi (t) , A 2 (t) , t G [0, T] be the rates of two Poisson inflows to the 
empty server, and 7 be the service rate. We want to estimate the difference 
b\ (t) — 62 (t) of the rates of the exit flows. We can couple the two service 
processes in the following way: let A (t) = min {Ai (t) , A2 (£)} . Then we write 
Xi (t) = X (t) + r/i (t) , where r/i (t) = A« (t) — A (t) . We will call the clients 
arriving with the rate A (t) as colorless, and we call the 771 (t) clients as red, 
while the r/ 2 (t) clients as blue. The colorless clients have priority in their 
service: if a colorless client arrives, then all the colored ones have to wait - 
even the one currently under the service. Then the difference \b± (t) — b 2 (t)\ 
is bounded from above by the sum of the exit rates of colored clients, which 
does not exceed 

I&! (t) - b 2 ( t )\ < 7 pr ( , ser ; er r is r?r d by a + , ) < 7 /Vi (o - a 2 (01 dt. 

1 1 w y colored client at the moment t J ~ ' J 1 v ; w 1 

Hence we have contraction with the contraction rate at most 7T, which is 
small for small T. We denote by X u (t) the unique fixed point of ty u . 

3. Finally we prove that the semigroup preserves the space of continuous 
functions. First of all we observe that the map depends on the initial 
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measure a in a continuous way. Therefore the same is true for X u (t) , the 
fixed point of ^/ u . Hence u (t) depends continuously on u (0) . ■ 

Let us consider the subspace C 2 C C° of functions /, which have the 
following properties: 

1. for every x G Z 5+ the function / has the first derivative ; 

2. for every x,y E Z 5+ the function / has the second derivative 9u r x \g u r y ) ] 

3. all these derivatives are bounded, uniformly in x,y. 

It is easy to see that the set of coordinate functions {u (x) , x E Z 5+ } on 
M. (K) can distinguish any two measures from M. (IK) . Due to the compact- 
ness of M. (K) we can apply the Stone- Weierstrass theorem, which implies 
that the subspace C 2 C C° is dense in C°. We now will show the following 

Proposition 3 For every t we have S t (Cq) C C 2 , where the subspace C 2 C C 2 
consists of all functions depending only on finitely many variables {u (x)} . 
In particular, the subspace C 2 is a core of the generator A. 

Proof. To do this we will use the Proposition 1 of the paper [DKVJ: 
Lemma 4 Consider the infinite system of equations 

j t z k (t) = a ^ (*) *i (*) + b k(t), t> 0. 

i 

Suppose that for all k 

^2\a ki (t)\ < a, \b k (t)\ < b exp{bt} , |^(0)|<c, 

i 

with a < b. Then 

\zk (t)\ < cexp {at} + - — — (exp {bt\ — exp {at}) . 

b — a 
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From (J7J) it follows immediately, that 
d f du t (v) 



dt \du (x) 



w£Y(v) 



du t (w) 
dun (x) 



duoM 



duo(x) 



with J2 we Y(v) K«< (*)l < a, 
so Lemma H] applies, and all the derivatives 
provided t < T. Further on, 



< 1, for some a < oo, uniformly in v, 
are uniformly bounded, 



du t (v) 
duo(x) 



d 



d 2 u t (v) 



dt \ du (x) du (y) 



^ a vw (t) f 



d 2 u t (w) 



w€Y(v) 



du (x) du (y) 



b v (t) 



with the same a vw (t)-s, while the term b v (t) , consisting of the products 
of the first derivatives ^fer ^fer , is also uniformly bounded, as was just 



duo(x) duo(y) 

shown, provided t < T. Therefore the derivatives 

bounded as well. 

For other functions we just use the chain rule. 



d 2 u t (v) 



du (x)du (y) 



are uniformly 



3.2 Equation for the evolution Vj^ and the conver- 
gence 

Now we will write the generator Am of the process W M . Let ip (■) be a func- 
tion on M.M- (In fact, we need it to be defined on a smaller set M.m H 
M. ^(-^Z 5 ) + j . Throughout this section the value of N will be fixed, and we 

will keep it just 1, in order to simplify the notation.) We will introduce the 
following notations for the increments of the measure v = (vq, v a , vg) (which 
are themselves (signed) measures on Z 1 or Z 2 ): 



^o,x 



1 if y = x 

(y) = { -i if V = x ~ 1 
otherwise 



x,y e Z 1 , 



&a,x a (y) 



1 if y = x A = {x A ,x BA ) 
-1 ify = (x A -l,x BA ) 
otherwise 



X4,?/ e z 
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1 if y = x A = (x A ,x BA ) 
&ba,x a (y) = { -1 if V = (x A , x B a - 1) , x A ,y e Z 2 , 
otherwise 

and similar definitions for the remaining measures A^ j- and A AB ^ X -. Then 



(AmV) (^) = ^c(^t;') [V («')-<0(v)] 



(9) 



xn>l 



'^ J± V A + A,3: m +(1,0> » U B ) _ ^ ( V 0> W A, Ufi) 



M ' 



M °, ^A, V S + B ' 3: f / +(1 ' 0) ) - V ^A, «b) 



+ Y M ^bav A (x a ) Y v o (x 



lp[v + A °a/° +1 , V A 



M ' 



v§) -il>(vo,v A ,v B ) 

+ Y M "Iabv b Ob) Y v o ^ ( v o + A °aT +1 ' v a,v§ - AA *f B J - if) (v , v A ,v s ) 

Xg.XAB>^ x O 

+ Y M ^av A (x a , 0) Y v b ( x b) *l> (yd, v A - AA f I Afi) , v B + Aab ' x s+(°^ \ _ ^ ^ v ^ v _j 

X^4>1 Xg 

Y m 1bV b {xb, 0) Y v a ( x a) if> (yo, V A + AbA 'm +(0,1) » v b ~ Ab 'm B ' 0) ) ~ *l> ( v 6, va, vs) 



XB>1 



Suppose now that the function ip is differentiable in each of the variables 
vq (xq) , %o £ ^\ V A ( x a) , X A e v b ( x b) ? x b Then each of the six 

increments in (|HD equals to the corresponding derivative ip'ofip, computed at 
some intermediate point. If moreover the function ip is continuously differ- 
entiable (which is implied by the twice differentiability), we can take a limit 
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as M — » oo, obtaining the convergence to the limiting operator 
(AiP) (v) 



= Y Oo) V A Oi) 

xo>l,x A 

+ Y : fvo(x )v s (xB) 

Xo>\,x B 

+ Y ^bava Oa) V£> Oo) 

X A -X BA >\,X 

+ Y "Iabv b Ob) Vq Oo) 

x B i x A :x A^l^ 

+ Y JbV b (x b ,0)v a (x a ) 

X A ,Xg.X B >l, 



dip(v) 



dip(v) 



+ 



dip(v) 



dip(v) 



vq{xq-1)) 8(vq(xo)) 8(v A (x A +l,x BA )) d(v A (x A ,x B A)) 



dif>(v) 



8ip(v) 



+ 



8ip(v) 



dip(v) 



"o(^O-l)) 8(vq(xo)) 8(v b (x b +1,x ab )) 8{v B (x B ,x AB )) 



8ip(v) 



dip(v) 



+ 



dip(v) 



{v {x o +l)) 8{vq{xo)) 8{v A (x A ,x BA -l)) 8(v A (x A ,x BA )) 
8ip(v) 8ip(v) , &i/i(v) dip(v) 



(i'o( x o+ 1 )) 9(vq(xo)) 8(y B (x B ,x AB -l)) 8(v B (x B ,x AB )) 
dip(v) dtp(v) , dil>(v) dip(v) 



8(v A (x A — 1,0)) 9(«^(a;A,0)) S(«g(i B «B+l)) S(»b(^b,ias)) 



8ip(v) 



8tp(v) 



+ 



8ip(v) 



8ip(v) 



v B (x B -l,0)) 8(v B (x B ,0)) d(v A (x A ,x BA +l)) 8(v a (x a ,x B a)) 



Let us apply this formula to the "coordinate" function ip (■) = (fi y (•) , where 
(j) y ( v ) = v (y) = v (yo) V A (Va) Vb (Vb) ■ The result is the one given by (jB} : 



(A(f) x ) (v) = -v (x 



v (x) 



( Y A «w) 

\a=0,A,B,AB,BA / 

Y la(l-S Xa )+-f A (l-S XA )5 XBA +-f B (l-5 XB )5 XAB ) 

\a=0,AB,BA J 



+ v(x- A ) (1 - 8 Xo ) (X AB (t) + Xba (*)) 

a=A,B 

+ v(x- A AB ) (1 - 6 XAB ) X A (t) +v(x- A BA ) (1 - S XBA ) X B (t) 
+ V (x + A a ) la 

a=0,AB,BA 

+ v (x + A A ) j A S XBA + v (x + A B ) <y B 8 x 



x AB ■ 



Since the space of functions of finitely many v-s coincides with that depending 
of finitely many u-s, that finishes the proof. 
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4 Fluid networks 



One of the key ingredients of the proof of out Main result is the investigation 
of the fluid (Euler) limits of various networks. They are introduced in the 
present Section. 

4.1 Fluid systems with one fluid 

The fluid systems with one fluid are the following simple dynamical systems. 
Consider the containers Vi, V 2 , V n , filled (partially) with water. Suppose 
that some pairs of these containers are connected by (directed) pipes, through 
which the water can flow. On every pipe i, j there is a pump p^- working, 
with the capacity of sending 7^ > units of water per unit time from Vi to 
Vj. The pumps are working constantly, and if the container Vi has less water 
than the pump p^- can handle, the result is that the pump sucks in whatever 
there is. For example, if the network is given by the graph 

Vx nv 2 n v 3 , 

with 7 12 = |, 723 = 1, then in a while the container V 2 will be empty, and 
the flow along the pipe 23 will be ~ (provided the water supply in V\ lasts). 
If a container Vi has several pipes ijk attached, then the water is shared by 
the pumps pij k proportionally to the capacities 7^. (this is relevant only in 
the situation when the level of water in Vi is zero, and all the incoming water 
immediately leaves it). 

Suppose now that at the moment T = all the containers are filled with 
water in the amounts of V{ (0) , and then we turn on all the pumps. We 
are looking on the levels V{ (t) , as they are changing in time. It turns out 
that there exists a time T' (depending on the system), after which the levels 
Vi (t) become stable and do not change anymore. (Some of them can in fact 
be zero.) In particular, such a system can never exhibit a cyclic behavior. 
Of course, the stability of the levels does not imply that the water in the 
network does not flow. It just means that for every container the amount of 
fluid entering is equal to the amount of fluid leaving. 

We will not provide the proof of this known statement (see [D] and the 
references there), since we will not use this fact. Below we will consider fluid 
networks which do exhibit cyclic behavior, and we find such examples among 
the fluid networks with several kinds of fluids. 
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4.2 Basic fluid network 



We will consider the following 5-dimensional dynamical system, A, which is 
a closed version of the open RS-network. It consists of three nodes: O, A 
and B, through which various fluids are passing. The node O has one type 
of the fluid, in the amount xo > 0, and that fluid flows into the nodes A 
and B in equal amounts. The rate of this flow 70 = 3, which means that 
three units of the fluid xo leave O per unit time (provided the supply lasts, 
of course), so each of the two nodes A and B gets § units of incoming fluids, 
A and B, per unit time. The amounts of these fluids are denoted by xa and 
Xg. The fluid A then goes to the node B, where it turns into the fluid AB, 
while the fluid B goes to A and turns there into BA. The corresponding 
rates are 7a = 7s = 10. The fluids AB and BA then go back to O, with 
the rates jab = Iba = 2. The last thing which has to be specified is the 
following priority: if the node A is in the state with both amounts xa and 
xba positive, then the fluid BA goes first. One can think that the fluid BA 
is heavier and of higher viscosity than A, so it goes to the bottom of the node 
A and flows out first (and relatively slow). The same applies to the node B. 




Basic fluid network. 



As stated, the system is not well-defined. For example, the dynamics is 
not specified if we try to start it from the configuration 

xa = a > 0, xb = b > 0, xba = xab = 0, xo = 1 — a — b. (10) 

The reason is that if the fluid A "starts first" , then it will create some amount 
of the heavy fluid AB in B, so the fluid B in B will be blocked. The same 
holds for the fluid B "starting first" . 
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The precise definition of the system, given below, follows [RStt \St2\ ID] . 
Consider first the simplest case of a single node with capacity 7 per unit 
time. Let the initial amount of fluid x(0) be given, an let Y(t) be the net 
amount of fluid that arrives to the node during the time interval [0,t]. In 
what follows we will call it the net inflow. The function Y(t) is monotone 
non- decreasing function, Y(0) = 0, and we assume that Y(t) is Lipschitz 
continuous, t > 0. It is easy to check that the evolution of the amount x(t) 
is given by 

x(t) = V(t) + U(t) } (11) 

where V(t) = x(0) + Y(t) — yt is called the virtual level of the fluid, while 
U{t) = max{0, — inf s£ [ 0ji ] V(s)} is the unused service capacity of our node. 
We introduce for the correspondence (fTTj) the notation 

x{-) = W{yx(0),Y{.)), (12) 

that is, W maps the initial fluid level and the inflow function to the evolution 
of fluid level. As we will show later (in Lemma [T3|) . the map W is Lipschitz 
continuous map from R x C[0, 00) to C[0, 00). 

Further on, let Z (t) be the total amount of fluid that leaves the node 
during the time interval [0, £]: Z(t) = x(0) + Y(t) —x(t). Again, the function 
Z (t) - the net outflow - is monotone non-decreasing Lipschitz continuous, 
with Z(0) = 0. By assumption, the derivatives z(t) = Z(t), y (t) = Y{t), 
existing a.e., satisfy 

= f 7 if x(t) > 0, 
\ V (t) otherwise. 

This property is the reason to call our discipline work-conserving; the server 
is always working at its full capacity. 

In the same way we can treat the node through which several fluids are 
passing. Thus we introduce the vector Y(t) = {Yi(t), . . . , Y n (t)} of the net 
inflows during the time interval [0,t], each Yi(t) being non- decreasing and 
Lipschitz continuous. The vector Z(t) will denote the corresponding col- 
lection of the net outflows. (Of course, it does depend on the priorities of 
the fluids.) Consider again the derivatives y%(t) = Yi(t) and Zi(t) = Z{(t) 
(they exist for almost all t and define Yi(t) and Zi(t) in a unique way once 
we fix Yi(0) and Zj(0) to be zero). Introduce also the workload rate by 
v(t) = YlVi^li 1 '1 where y are the service rates. The service discipline of 
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our node is work- conserving, if the following property holds: once ||x(t)|| > 0, 
we have 

5>(*hr 1 = i; (13) 

otherwise 

z>(t)=yi(t) (14) 
for all i. The following statement is immediate: 

Proposition 5 Let x(0) = and v(t) < 1 for almost all t > 0. Then 
x(t) = 0, t > 0. 

The system described in the beginning of the present subsection corre- 
sponds to the following specification of the above general formulation. We 
have Y(t) e (IR 5 ) + , and in terms of the map W (see (fl2l ) the evolution is 
given by the equations: 

x o (-) = W(3,x o (0),Yo(-)), (15) 

xba(-) = W(2, x BA (0), Y BA {-)), (16) 

x A {-) = W^(10, hx BA (Q) + x A (0), hY BA {-) + Y A {-)) - 5x BA (-), (17) 

the equations for x B and x AB being symmetric. Then for the derivatives, 
whenever they exist, we derive from f|T5l) - f|T7j) that 



d , , 
dt Xo (t) = 

JXAB (t) 



dt 



x A (t) 



Uo (t) — 3 if xo (t) > 0, or if x (t) = and yo (t) — 3 > 0, 
if x (t) = and y (t) - 3 < 0; 

(t) - 2 if xab (*) > 0, or if x AB (t) = and y AB (t) - 2 > 0, 
if x AB (t) = and y AB (t)-2< 0; 

if x BA (t) > 0, 

if xa (t) > 0, x BA (t) = 0, y AS (t) - 2 < 0; or 
if xa (t) = x B A (t) = 0, 2/ab (t) - 2 < 



VA (t) 

VA (t) - io^w 



and ?/a (0 — 10 



2-y AB (t) 



> 



if xa (t) = s B A (t) = 0, y AB (t) - 2 < 



and ?/a (£) — 10 



< 0: 
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the equations for xb and xba being symmetric. 

In fact, for most points in the phase space the above differential equations 
are sufficient to describe our dynamical system, since the set of values t 
where one of the derivatives does not exist is nowhere dense. However, this 
is not true for all points, and so we need to use more complicated equations 

(Ha-LTZD. 

The above relations define what we will call Non- Homogeneous Dynamical 
System, NHDS. We will denote this dynamical system by A (Y) , since it 
is driven by the inflow Y. This is just a usual non-autonomous dynamical 
system. All the non-linear dynamical systems, which will appear below, 
correspond to different choices of the flows Y. 

Let now x (t) be the trajectory of the NHDS, corresponding to the initial 
state x (0) and the given inflows Y (t) . We define the closed fluid network 
evolution A of the point x (0) as the evolution x (t) under any dynamics 
A [Y (i)J , for which the following relations between the inflows Y (i) and 
the outflows 

Z(t) = Z^{t) = x(0) + Y (t) - x(t) (18) 

hold: 

Yo (t) = Z% it) + 2%\ it) , (19) 
Ya (t) = \zf (t) , (20) 

Y AB (t) = Zf (t) , (21) 

and symmetric relations for A and BA variables. The set of all solutions 
Y (t) of the equations will be denoted by y (x (0)) . 

It is well known that the set y (x (0)) is not empty for any initial state 
x(0). We reproduce here the proof, since analogous argument will be used 
throughout the paper. 

Proposition 6 For any point x(0) there exists at least one trajectory of 
closed fluid network evolution A, passing through it. 

Proof. It suffices to prove that a solution exists in any given bounded time 
interval [0, T\. For a given x = x(0), consider the map G : Y(-) — > Z X,Y (-). It 
is clearly a continuous map of C[0, T] into itself. The outflow Z(-) is Lipschitz 
continuous (by Lemma [T51 below) . with Lipschitz constant L independent of 
x(0) or ?(■). Hence G takes the convex compact set of L-Lipschitz continuous 
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functions Y(-) on [0, T], satisfying Y(0) = 0, into itself. Therefore, by the 
Brouwer theorem, the map G has at least one fixed point. ■ 
The system A has the following properties: 

• The total amount of fluid | x (t) \ = xo (t)+XA (t)+...+XB (t) is evidently 
conserved. We will assume \x (t)\ = 1. 

• For some initial states x (0) the equations (IT31) - (|2"T1) have multiple so- 
lutions. This is true for all initial conditions x (0) , given by (1101) . The 
equations (fl9|) — f l2T|) have in this case three solutions, so there are 
at least three dynamical systems, defined by x (0) . The first solution 
corresponds to the flow rates y A B = 10, Uba — 0, the other one is 
symmetric to the first one: y B A = 10, uab = 0, while the third one 
is given by the flow rates AB -> O, BA — > 0,A -> AB, B -> BA all 
equal to 1?, O — > A, O — > 5 equal to |. One can say, having that 
property in mind, that at some points the uniqueness of the trajectory 
breaks down, so it can happen that for two trajectories x' (t) , x" (t) 
we have x' (t) = x" (t) for t < t , but x' (t) ^ x" (t) for t > t . That 
just means that there are two different (non-autonomous) dynamical 
systems, which have trajectories, coinciding for t < t , but not for 
t > t . 

• The curve x (t) = {1,0,0,0,0} is a trajectory, i.e. the point * = 
{1,0,0,0,0} is a fixed point of A. Its flow rates are constant; their 
values are: Ua = Ub = Vab = Uba — §, Do — 3- Let us check that with 
these flow rates the amount of the fluids in A and B will stay zero. 
Indeed, — + = ^ + ^ = ^<l. So our claim follows from the 

' 7A "IBA 10 2 10 

Proposition [5j In fact, the nodes A and B are even underloaded, which 
means that in the long run each node gets on the average less fluid than 
its serving capacity is. Note also that there are trajectories x (t) such 
that x (i) = * for t < t , but x (t) ^ * for t > t , with any t . 

• There are (not necessarily uniqueness) points in (IR 5 ) + , from where 
(some) trajectory goes to *. One such family is the set of points U = 
{x : < xa = xb < \i xab = %ba = 0, xo = 1 — xa — xb} , and the flows 
are: AB -> O, BA -> O, A -> AB, B ^ BA all equal to f , O -> A, 

O — >• B equal to |. In fact, from every point in U two more trajectories 
start. The values of the flow rates for one of them for small initial seg- 
ment of time is given by: y A = y B = |, y AB = 10, y BA = 0, y Q = 2. It is 
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describing the situation when the light fluid A flows to the node B, and 
the resulting heavy fluid AB in the node B blocks the fluid B in the 
node B from exiting. The second solution is obtained by interchanging 
A and B. There is a bigger set U D U, having dimension 2, starting from 
where one can get to *, making on the way some choice of the trajectory; 
U = {x : x A = x B , x AB = x BA , x = 1 - x A - x B - x AB - x BA } . 

• There is a cycle C C (M 5 ) + , such that if x (0) G C, then x (t) G C for all 
t > 0. For example, the point {x : x A = 1, x B = x AB = x BA = Xo = 0} 
belongs to it. All the points of the cycle C are uniqueness points. By 
Tc we will denote the time to go around the cycle C once. We will now 
describe, just applying the definitions in a straightforward way, the 
cycle C, started from x (0) = {x A — l,x B — x AB = x BA = xo = 0} . 
The first part of it happens for t G [0, |] , during which time the 
component x A (t) decays linearly from the value 1 to 0. The component 
x B (t) grows linearly, and x B (|) = |, the component x AB (t) grows 
linearly, and x AB (|) = |, while x BA and xo stay zero. On the next 
segment, t G [|, l] , the component x AB (t) decays linearly, with rate 
1, and so it vanishes at t — 1. The component x B (t) continues to grow 
linearly, with the same speed, so x B (1) = 1. Other three components 
remain empty. So at time 1 we find ourselves accomplishing one half 
of the cycle. Therefore Tc = 2. 

• If x (0) ^ U, then x (t) G C once t>T, where T = T (70, 7a, 1b, Iab, 7b a) ■ 
We will prove this claim under additional assumption that dist(x (0) , C) = 
inf^ec J2i \ x i ~ Ui\ < e f° r some small e, since this will be sufficient for 
our purposes. 

Lemma 7 Local attractor. There exists an e > such that for any ini- 
tial point x (0) with dist(x (0) ,C) < e we have x(t) G C once t > T = 

T ' ino,lA,lB,lAB,lBA) ■ 

Proof. Since for every point c on C we have either c AB = or c BA = 0, we 
see that at least one of the coordinates x AB (0) or x BA (0) has to be less than 
e. Let us start with the case that both of them are positive, and we can then 
assume that x AB (0) < e, x BA (0) > x AB (0) . Then for a short initial segment 
of time, [0,ti] , both coordinates x AB (0) and x BA (0) decay with the same 
constant rate 2, until x AB will vanish, which will happen at the moment 
ti = 2A|M < 1 Also, for every point c on C we have cq = 0. Therefore 
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x (0) < e, hence x (£1) < e + | = y , since the total arrival rate to O does 
not exceed 4, while its service rate equals 3. 

Now consider the case when xba (^i) > 0. While xba{^) keeps being 
positive, there is no flow from A to B and, as a result, no flow from B to O. 
So if xba if) stays positive for £ e [£i, £2 = ^i + (*i)] , then i G (£) becomes 
at £ = £2 ~ since it decays with rate 1 when xab = and xba > - and will 
stay while xba (£) is positive. The coordinate xab (£) stays zero as well. 
When finally xba (t) vanishes for the first time, at £3 > t 2 , the value xb (£3) 
has to be already zero. Which means that xa (£3) = 1, so x (£ 3 ) G C. 

In the remaining case, the first point, £3, where the coordinate xba (t) 
vanishes, belongs to the segment [t 1: t 2 = £1 + x (£1)] • (That case contains 
the situation when xba (ti) — 0.) Since the time £3 < 2e, and since initially 
the point x (0) was close to C, the same is true for the point x (£3) . Since 
Xab (h) = xba (h) = 0, we conclude that one of the two coordinates - either 
Xa (£3) , or x B (£3) - should be e-close to 1, while the remaining one, as well as 
%o (£3) , should be e-close to 0. Suppose xa (£3) ~ 1- Note that the evolution 
of the point x (£3) is not uniquely defined if both ^(£3) and xb (£3) are 
positive. As was explained above, there are, in fact, three options to choose 
from: 

i) If the choice is that the fluid A "goes first", then after a small time the 
coordinate xo will vanish, and after the time of order | the coordinate xa 
will vanish as well, and we find ourselves on C. 

ii) If the fluid B "goes first", then after a small time (of order e) first xb, 
and then xba will vanish, and the fluid xa - which is in the amount of order 
1, will start to decay, so we find ourselves in the situation just considered. 

Hi) The remaining option is when both fluids A and B "go simultane- 
ously" into, resp., AB and BA, at the rate As was explained above, the 
result will be that the levels of the fluids A and B will decay with the rate 
— I = |, while the level xo will be correspondingly raising. But no later 
than the time 6e the level xb will vanish, and since xab and xba both were 
already zero, we again are in the situation considered above, with xo being 
of order e. In fact, at any time when the system is in phase Hi), it has three 
options: to pass (forever!) to the phase i) or ii), or to stay in phase Hi). The 
above arguments in i) — Hi) stays valid in this case and we get the required 
assertion. ■ 

Warning. Our analysis shows that the time T needed to reach the cycle 
does not vanish with e. 
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Continuity. The above proof implies the following (weaker) substitute 
for the property of the continuous dependence of the trajectory on 
the initial condition. Let x (0) and c (0) be two initial points, and 
suppose that |x (0) — c(0)| < e, and c(0) G C. Then there exists a 
constant C, such that for any t and for any version of the x-trajectory 
we have \x(t) — c(t)\ < Ce. The condition c(0) G C is, evidently, 
crucial; without it our statement fails. 

The cycle C depends on 7o, 7a, 7b, 7ab, 7ba, and is non-trivial for our 
choice of these parameters. For some other values of 70, 7a, 1b-, Iab, Iba 
it is reduced to the point *, which then is a stable fixed point. 

All the above properties of our system would still be valid if we perturb 
slightly the vector 7 = {70, 7a, 7b, 7ab, Iba} of the parameters around 
the point {3, 10, 10, 2, 2} of our choice (even if the perturbation does 
not respect the symmetry ^ A = ~/ B , Iab = Iba)- 

Let x (t) C C be a cyclic trajectory. Let us denote by A c (t) the 
corresponding (periodic) function of the inflows. The stability prop- 
erty just formulated implies immediately that our system in the cyclic 
regime is underloaded. In other words, for {70, 7a, 7b, Iab, Iba} = 
{3, 10, 10, 2, 2} there exists 5 > such that we have 



1 <- T ' 



lo Jo 



° \ C o(t)dt< (1-5)T C , (22) 



— / * \ c A (t)dt+— [ C X c BA (t)dt<(l-5)T c , (23) 

1A Jo IBA Jo 

and the same relation for the B node. 

We now want to consider the "open" system A°, which is obtained from 
A by the following construction: every exiting fluid now goes not to the 
corresponding node, but leaves the system. On the other hand, there is some 
inflow, A (t), entering the system from the outside. 

For the future use we introduce now the following compact subset /C = 
K (7) C (M 5 ) + : 

K = \xe (1R 5 ) + : max < -xab H %b, -%ba H xa, —xo V < 10 > . 

\ K ' \ 2 10 2 10 3 J J 
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Lemma 8 Let the open system A° be in the state x (0) G (IR 5 ) + . Consider 
the functional L onM 5 , given by 

L(x) = \ ;1 x \ ii, ^ G/C ' (24) 

L max | ^ab + jqX B , 2 x ba + jqXa, 3^0) otherwise. 

Suppose that L(x(0)) > 10. Then there exists a constant C > 0, such 
that for all external flow rates {X(t) ,t G [0, Tc]} which are close enough to 
{X c (t) ,t G [0, Tc]} in the L 1 distance, we have 

L(x (0))-L(x (T c )) >C 

Proof. For the case of the inflows with rates A c (t) our statement follows 
from the underload property, due to the relations fl22l) - fl23l) . Therefore it 
holds for inflows that are close enough, by continuity. ■ 



4.3 M coupled fluid networks 

Let Am be the dynamical system on (R 5M ) + , obtained from M copies of 
A, interconnected in the mean- field manner, as follows. Each node Oi, 
i = 1,...,M, is connected to all of the nodes Aj, Bj, j = 1,...,M, and its 
fluid, of the amount xo,i, flows into the nodes Aj, Bj in equal amounts. The 
rate of each of these flows is now M _1 ^y- = ^f^, which means, as before, 
that three units of the fluid xo,i leave Oi per unit time, so each of the set 
{Aj} and {-Bj} gets | units of incoming fluids, A and B, per unit time. In 
a similar way, the fluid A from every node Ai is splitted among all nodes 
{Bj} , so the rate of every individual flow A4 — > Bj is M~ lr jA = , and so 
on. The priorities are kept the same: if the node A iy say, is in the state with 
both amounts xa,i and xba,i positive, then the fluid BA goes first. 

Again, we first describe the open network, i.e. the Non-Homogeneous 
Dynamical System. We will not need the general case here; it is enough for 
us to consider the net inflow defined by the same function Y(t) G (1R 5 ) + as in 
the previous subsection; every node then gets jj-th part of the inflow, so, for 
example, for each i = 1, M the net inflow function of the node 0{ equals 
to jjYq (t) . Once we are also given the initial values x G (M 5M ) + of the fluid 

levels, the net outflows Z X J* it) , a G {O, A, B, AB, BA} , i = 1,...,M, are 
defined as above, see (1181) . 
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Passing to the closed system, instead of relations (IT5|) - (12T|) . we impose 
the relations 

M 

Yo(t) = Y,( Z Aii(t) + Z*Blm (25) 
1=1 

- 1 - 

y *j (*) = E 2 Z % {t) ' (26) 

1=1 

M 

Yabj (*) = £ Z l? (*) • ( 27 ) 

i=l 

Again, all the functions and Z#(t) are Lipschitz continuous and, hence, 
differentiable almost everywhere. These derivatives will be denoted, again, 
by y*(t) and z*(t). For each node (say, Ai) of Am, the evolution of the amount 
of fluids XA,i(t) and XBA,iif) is found from the corresponding inflows -jgY^(£) 
and jjYsAit) and initial states XA,i(0) and xba,i(0) in the same manner as 
for the network A. 

We will consider trajectories x M (t) G (R 5A/ ) + with \x M {t) \ = M, i.e. we 
have the unit amount of fluid per elementary system A C Am- 

For every M we will define now another dynamical system, acting on 
(R 5 ) + . This one, also denoted by Am, will be of central importance for 
the present paper. However, it will be not the usual dynamical system. 
It will be defined not as a group of transformations of (R 5 ) + , but directly 
on the (sub)set Mm of some atomic probability measures on (R 5 ) + . This 
transformation will not be linear on M, and for that reason we will call it 
non-linear dynamical system. In fact, it is just a convenient representation 
of our initial dynamical system on (R 5M ) + , suitable for passing to the limit 
M — > oo. The construction is very simple: 

To every point x M = { (x Q , x\ = (x A , x l BA ) , x* s = (x l B , x l AB )) , % = 1, M} G 
we can assign a probability measure on (R 5 ) in the following way: 
we put no = jg E]=i l*A = Ji Ei=i 5 4> f^B = if E]=i 8 x%, and we de- 
fine \i = fis M — Ho x Ha x Hb e M-m- Now, if the point x M evolves according 
to A M , so is the measure \i\ moreover, the evolution fi (t) of fi is well defined 
and does not depend on the choice of the preimage, so if \i.m x = /z_m 2 , then 

x l x 2 

v-xf 1 (*) = ^4' 2 (*) • 

The set of measures Mm on (R 5 ) + consists of all measures /z, having the 
properties 
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1. /x is a product, 

fi = (p Q , pa, a*b) E ^x/iix/j 5 En 5 [/i] x n A [/i] x n 5 [p] , (28) 

of probability measures on M 1 = {xo} , resp. M. 2 = {xa,xba} and 
= {^b,^ab} • Here we denote by H^-s the various projections (or 
marginals) , 

2. we have J \x\ dp = f xodpo + / (xa + x ba) dp a + / ( x b + £ab) dp§ = 
1, 

3. we have /x = tf X/i^i 4, f° r some (not necessarily distinct) ^ G M 1 , 
i = 1, M, likewise Pa=m~ E^i &b<> = E 2 =i ^, G 

Properties of Am : 

• The set of the fixed points of consists of measures /x = U d [p] x 
5^=0 x 5 Xg =o, where <5 XA =o and 5 XB= o are unit atoms at the origin, while 
Ho [p] is the projection on the coordinate Xo- In words, that means that 
all the fluid stays permanently in the O-nodes (in arbitrary amounts, 
adding up to M). This fact follows from Proposition [51 

• If p (0) = 5 S G M-m, then /x (t) = S S (t), where x (t) is the trajectory of 
A with x (0) = x. In particular if x G C, then x (t) G C. 

• Note that the dynamics Am on A4m is "non-linear" , in the sense that in 

general in the situation when /x (0) = a/x' (0) + (l — a) /x" (0) , /x (0) , pi (0) , /x" (0) G 
M M , < a < 1, we have /x (i) 7^ a/x' (i) + (1 - a) p" (t) for t > 0. 
(There is nothing strange or unusual in this relation, since the dynam- 
ical system in question is itself defined on the space of measures as its 
state space, and not on the (M 5 ) + .) 

• Let Pkrov be the Kantorovich-Rubinstein-Ornstein-Vaserstein distance 
on the probability measures on (M 5 ) + , corresponding to the metric 
p(x,y) = Yli=i\ x i~Vi\ 011 (^ 5 ) + • (We recall briefly, that if /x,/x'are 
two probability measures on a metric space (X, p), then Pkrov (p, p') = 
inf K J p (x, x') dn (x, x') , where the inf is taken over all probability mea- 
sures k on X x X, such that k(A x X) = p (A) , k(X x A) = p' (A) .) 
Suppose that the initial measure p (0) is close enough to the cycle C, 
which means that for some x G C we have 

Pkrov {p (0) , 6 X ) < e. (29) 
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Let Y (t) G y (fi (0)) be one of the possible net inflows, corresponding to 
the initial state \i (0) , and \i (t) be the corresponding evolution. Then 
there exists the time T = T (M, 70, 7a, 1b, Jab, Iba) , such that for all 
t > T we have jj, (t) = S^a) with x (t) G C. Moreover, let us define the 
compact Ax by 



A K = {x = (x Q , x A , ...,x B ) : x < K, x A < K, x B < K} 



and let \x 



K 



(t) be the evolution of the restriction \x (0) 



under 



K 



the same evolution, defined by the flow rates Y (t) G y (fi (0)) . (Once 
the inflows Y (t) are fixed, the evolution becomes the usual (non- 
autonomous) linear dynamical system, so we can apply the dynamics 

to the summand /1 (0) of the measure // (0) .) We choose K to be 

large enough, so that from Pkrov (z 1 (0) ,8 X ) < e for some x G C it 
follows that 

fi (0) [K]>1- e. (30) 
We now claim the following: 



Lemma 9 Under conditions <\29t\ and (jffflj) , there exists the time moment 
T' = T' (K, 70, 1a, 1b, Jab, Iba) , such that for every t > T' (uniformly in 
Ml) there exists a point x (t) , such that \i it) is just the atom at that 



point: fi 
while c - 



K 

(t) = c5$u)- Moreover, dist(x (t) ,C) < c, and c 
as K — > 00. 



ii (0) [K] 



Notes. 



• We will prove only the statement about the existence of the time mo- 
ment T', since below we will not use the time moment T (M, 70, Ja, 1b, Jab, Iba) ■ 

• Our proof can be extended literally to the case M = 00 of the next 
Subsection. 



Proof. Note first that due to the mean-field nature of our graph, the flows to 
all the A-nodes are equal at every time moment (as well as to all the .B-nodes 
or O-nodes). Consider now any subset Q of the A-nodes, \Q\ < M, and let 
Iba be the index, for which (xba)j ba (0) > (xs J 4) i (0) for all i G Q. Then, 
clearly, this relation holds at later moments, i.e. (xba)i ba (t) > ( x ba)i (t) ■ 
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In the same way, define the index I a as the one, for which (xa) Ia (0) > 
{xa\ (0) . Then, if all the variables {xba)i (0) are equal for % G Q, the relation 
(xa)j a (t) > (xA)i (t) holds at all later moments. We will use this property 
for the set Q = {i} of indices, which satisfy (xba)i (0) < K, (xa)^ (0) < K. 

Let us show that there exists the time moment, such that before it every 
node in Q will be empty for some time duration. In view of what was said 
before, it means that after that time moment all the nodes in Q will be 
synchronized. To see this depletion, note that the initial supply of the fluid 
BA at any node i G Q is not exceeding K, so it will be over before t' = -j. 
Next, there exists a moment t", after which ^ (say) of "atoms" of the heavy 
fluid BA, passing through our node i G Q were at earlier moments at some 
O-node. Indeed, another option would be that such an atom was staying 
at some B node for all the time duration t". That means that the initial 
total amount of fluid at this node was very high, once t" is chosen to be 
large. However, the proportion of such nodes has to be small, due to the 
simple fact that the total amount of fluid per node is of the order of one. 
But the rate, at which the fluid goes from the O node to the B node is never 
higher than |. Let K" be the amount of fluid at our node i at the moment 
t" . Clearly it is at most 12t" . Suppose the node is not empty during the 
time interval [t",t" + T] . Than it works all the time at full capacity. Let 
< k (t) < 1 be the fraction at moment t of the capacity of the node, used 
by the heavy fluid, while the remaining fraction 1 — k (t) is used by the light 
fluid. Then the amount of heavy fluid, which left the server during this time 
interval, is 2 J „ k (t) dt, while the corresponding amount of the light fluid 

is 10 J t „ (1 — k (t)) dt. Since the light fluid flows into the node with the 
rate at most |, we have that the relation 

f t"+T o 

10 / (l-k(t))dt< -T + K" 
Jt" 2 

has to hold, since the amount of light fluid, leaving the node, can not exceed 

the initial amount present at the node plus the amount which came to the 

node during the time interval T. For the heavy fluid we similarly have 

2^ (t) *<fr + io( I Lr) + A-». 

The two relations imply that 

10T < ( 9- ) T + 6K". 
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So T < 12K", which establish our depletion claim for the A (as well as for 
B) nodes at some moment t'", independent of M and e, provided only that 
e is small. 

Thus far we were not using the condition pkrov (/-£ (0) , S x ) < e, without 
which our claim about the existence of the time moment T' (K, 70, Ja, 1b, Iab, Iba) 
is not valid - see, for example, the first property of the dynamics A^. We will 
use it now, in dealing with the O nodes. Due to the above discussion and the 
continuity property, our statement is reduced to the following one: consider 
the initial measure \i (0) , having the properties: Pkrov (/•* (0) , Sx) < s'{= Ce, 
see Continuity Property of the previous section), while in the decomposi- 
tion \i (0) = n (0)0 x n (0)^ x //(0)j we have \i (0)^ = {\ — e)S XA + xa, 
/i(0)g = (1 — e) S$ g + Xg, with the vectors xa,x b ^ 2 close to the corre- 
sponding projections (xa,%ba), resp. (xb,xab) of the vector x. But that 
means that the flows into the O nodes will be almost always almost equal to 
these on the cycle C, so in finite time all of the O-nodes which initially have 
their levels < K will become empty. Indeed, on the cycle the flow to the O 
node has rate 2, while the capacity of these nodes equals 3. 

Summarizing, we have thus far that after a finite time, independent of 
M, all the nodes in Q are synchronized, and moreover all the O-nodes in Q 
are empty. The application of the Continuity Property and the Attraction 
Lemma [7] finishes the proof. ■ 



4.4 M — > 00 fluid network 

This is again a dynamical system, Aoo, acting on probability measures M. = 
J\A ^(M 5 ) + j on (M 5 ) + . One way of defining it is to say that the family // (t) G 

M. is a trajectory of A^, iff for any M there is a trajectory /i M (t) G M.m 
of Am, so that for every t we have \iu (t) — > n if) weakly. 

Now we will give another description of A M , which does not make use of 
the limit M — > 00. As was the case with the dynamics A and Am, we will 

define first the NHDS version of A^. Let the measure /io G M. ^(M 5 ) + j and 
let the function Y (t) G (1R 5 ) + be given, which is the net inflow of our fluids. 
Then the corresponding evolution // = G M. ^(1R 5 ) + ^| of the state /i 
is defined to be just the evolution of the measure /xq under the dynamics 
A (Y) . Again, the net outflows (t) , x G (JR 5 ) + are given by ((THD . 

Now we can define the evolution [L — G M. ^(M 5 ) + ^ | of the initial 
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measure /j,q under Non-Linear Dynamical System A^ (NLDS) as the NHDS 
evolution of it under any dynamics A (F) , with Y satisfying the equations: 

Y (t) = f{zf B {t) + zf A {t))dvo{x), 

Y A {t) = \ j zf{t)d^{x), (31) 

Y AB it) = j Zf (t) dfio (x) , 

and symmetric relations for B and BA variables. The set of all such flows 
Y (t) will be denoted by y (/x ) • 

We are calling the dynamical system Aoo non-linear, since in general we 
will have for /x = \ (^'q + h'q) that ji t ^ \ (n' t + Z-4') when t > 0. Note also 
that if /i = <5 5 (o) for some point x (0) , then the family ji t is a A^ trajectory 
iff fit = Sx(t), with x (t) being some A-evolution of x (0) . 

Properties of A^ : 

1. the set of the fixed points of A^ consists of measures /j, = U 5 [/j] x 
5 Xa=0 x S x - =0 , where U 5 \p] is the projection on the coordinate x<> 

2. if /i (0) G Mm for some M, then the A M -dynamics and A^-dynamics 
with that initial data coincide. In particular, if ji (0) = 5 X G M. with 
x G C, then \x (t) = 5 x ( t ), where x (t) is the trajectory of A with x (0) = 
x. 

For all our purposes it is sufficient to prove the following property of our 
NLDS. 

Proposition 10 Let the measure fi = // (0) on (IR 5 ) + have the following 
properties: 

i) Unit mass: 

/ (xa + x B + x B A + x A b + x ) d/i = 1, (32) 

J(R 5 )+ 

ii) Exponential moment condition: for some a > 0, A < oo we have 

(exp{aL(x)}) m <3A (33) 
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(see 

in) For some x G C we have 



Pkrov O (0) , S x ) < e, (34) 

with e small enough (depending on a and A). 

Consider now some Non-Linear Dynamical System (NLDS), defined 
by the initial state fi (0) . In other words, A^ = Aoo (Y (■)) > f or some Y (•) G 
y (/i (0)) . Then for t — > oo iae evolving measure \x (t) satisfies : 

I) 

Pkrov (p (*) , -> (35) 

for appropriate z(t) G C , 
II) 

(exp{aL(x)}) MW ^0. (36) 

Moreover, the convergence in f lffjp , f lffffi) zs uniform over all Y (■) £ y (fi (0)) 
and all initial measures [i (0) satisfying . 

To establish it we first prove the following simpler fact. 

Proposition 11 For every T,e there exists a value e(T,e) , such that the 
following holds: 

i) For every T we have e (T, e) — > as e — > 0. 

m) Lei /i &e any measure satisfying i) - iii) above. Then for every t G 
[0, T] we aawe p^ov (a* (*) > ^(t)) < e ( T , e) • 

Proof. Let the sequence of measures \i n (0) converge to 5 X in the KROV 
metrics, and let it satisfy the conditions of the Proposition [TU1 Let us consider 
the set of trajectories (t) , < t < T. The family of these trajectories, 
viewed as functions of t G [0, T] , is a family of uniformly bounded functions 
(due to the compactness of the set of measures with the properties i) - Hi)), 
which also are equicontinuous. Indeed, every version of the vector field, along 
which any of the measures \i n (t) has to evolve, is continuous and bounded 
in norm by the constant 7^ = 10. Therefore it is compact, and so has a limit 
point, the function \i (t) , with fi (0) = 5 X . But fi(t) has to be a trajectory 
of NLDS, and since there is just one such trajectory starting from 5 X , we 
conclude that \i (t) = 5 x t t y 
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Therefore p n (t) — > 5 x ( t ) for every t G [0, T] , and this convergence is 
uniform in t. The existence of the function e (T, e) follows from the compact- 
ness of the balls {p : Pkrov (p, $x) < a } of measures satisfying the moment 
condition. 

Proof of the Proposition 1101 

1. To begin with, we prove that once e is small enough, there exists the 
time moment Ti, at which for any initial p (0) , satisfying conditions of our 
Proposition, 

(exp{aL(x)}) MTi) < e. (37) 

Indeed, suppose first that the value of the exponential moment of the initial 
measure p (0) is fixed, (exp {aL (x)})^ = M, and also that 

Pkrov {p if) > ^z{t)) < £ for all t G [0, T c ] , with e small enough, (38) 

where Tq is the time it takes to go once around the cycle C. We claim that 
at the moment T c we have (exp {aL (x)})-^ < cM, where c < 1 is some 
constant, which depends only on the parameters of our model. In particular, 
c is independent of jl (0) . 

To see that we note first that for any point x with L (x ) > 10 (see 
definition (jMD) and under Cyclic Dynamics = (Y (■)) , for Y (•) e 
y (S z ) , z G C, we have after the time shift by T c that L (xt c ) < L (^o) — 
C (7*) , where the constant C (7*) > depends only on the service rates 7*. 
This follows directly from the Lemma [8] of the Section 14.21 From the same 
Lemma and due to the condition (1551) we have that under any dynamics 
Aqo (Y (•)) with Y (■) G y (p, (0)) we have for the same point that L (xt c ) < 
L (xq) — \C (7*) . So we have proven our claim, with c = exp { — |C (7*)} . 

Let now k be the smallest integer, such that 3Ac k < e. We want to repeat 
k times the procedure of the previous paragraph. Its duration is, evidently, 
kTc- Suppose that the parameter e is so small that the function e(kTc,e) , 
defined in the Proposition fTTl satisfies e(kTc,e) < e (see ( 1551) ). Then the 
desired repetition is possible, and so (1571) indeed holds, with 7\ = kTc- 

2. Summarizing, at the moment T\ we have: 



(exp{aL(x)}) M(Ti) < e, (39) 
H) Pkrov {p (Ti) , 5 Z ) < e (T 1} e) for some z eC. 

Let us use now the compact K. From ( 1591) immediately follows that 
p (Ti) [JC] > 1 — e. Hence, due to the Lemma[9]of Section fl~51 there is the time 
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moment T 2 = T 2 (/C, e) , after which the restricted measure p (Ti) would 

evolve under the dynamics Aoo (Y (■)) with Y (■) E y (p (0)) to a 5-atom, of 
mass at least I — e. According to Lemma [TT1 - or rather its M = 00 version 
- for some z G C the distance Pkrov (p (Ti + T 2 ) , £ 2 ) at that moment would 
be at most e (T 2 , e (Ti, e)) , which is small. If it would have been the case 
that this atom has unit mass, than after time Tq it would already be on the 
cycle C. Since, however, we know only that its mass is above 1 — e, we can 
claim that at the time moment Ti + T 2 + Tq its distance from C is at most 
Kie, where K\ is some universal constant. Therefore for the total measure 
p (Ti + T 2 + Tq) we can claim that 

ii') Pkrov (p (Ti + T 2 + T c ) , <5 2 ) < K 2 e for some z 6 C and some univer- 
sal K 2 , while 

i') (exp {aT (^)}} m (t 1 +t 2 +t c ) < ce ^ w ith c < 1 the same as in the first step 
of the proof. 

So what happened is that the estimator e of the exponential moment, 
appearing in (j3U|) , multiplied by K 2 , is the estimator of the KROV distance 
on the next step, while the estimate of the exponential moment is improved 
by a constant c. This argument can be iterated; on the next step we will have 

i") (exp {aL (x)}) ^ {Ti+2(T2+Tc)) <c 2 e, 

H") Pkrov (p (Ti + 2 (T 2 + T c )) , 5 2 ) < i^ 2 c£: for some z e C, 
and so on, which completes the proof. ■ 



5 Euler scaling limit of the NLMP 

In this Section we will show that in the limit of high load our Non-Linear 
Markov process tends to the fluid model, introduced above. 

Theorem 12 Suppose that for every N we are given the initial state v 
of the (scaled - see (fTJ)^ NLMP V^, and the sequence v converges to the 
measure p in the KROV metric, i.e. Pkrov (v n , p) — > 0. Consider all the 
KROV -limit points p (t) of the set of trajectories {P N (t) = u N (Nt) , t G [0, T]} , 
iV = 1, 2, ... . Every such limit point, called Euler fluid limit of the NLMP, 
is necessarily a trajectory of the fluid model A^. 

(Any trajectory of A^, obtained via this limit, will be called a fluid solu- 
tion.) 

Note. As was mentioned above, the system A^ does not possess the 
uniqueness property. The uniqueness property for the subclass of the trajec- 
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tories of Aoo - the fluid solution trajectories - does not hold as well. That 
means that the above set of all limit trajectories /i (t) might contain more 
than one element, and so the trajectory ji u N (t) , which satisfies h u n (0) = /x, 
does depend on the sequence u N — > /i. 

Proof. 1. We start by considering the open system. Then we can 
consider every node separately. Let us take the node A, say. In the open 
system case its evolution is defined by prescribing the initial state, - the 
distribution of the quantity q\- (0) G 1R 2 , which is the initial queue, scaled by 
the factor - together with the rate function A^ (t) = {A^ (t) , A^ A (t)} G 
M 2 , t > 0, which defines the Poisson flows of incoming clients. For our 
applications it is enough to consider the case when all our rate functions A-s 
are uniformly bounded: 

Af (*) < CN, (40) 

since this is definitely the case for our closed system. The service times of 
the clients are exponential, with respective rates N'-fA, N'-fBA (this scaling 
is due to the Euler limit we are going to study). As above, the BA clients 
have priority. That means that if, while a user of class A is being served, a 
class BA user arrives, the service of A user is interrupted until the moment 
when there will be no BA users in queue (the preemptive priority service 
discipline). 

For the future use we will introduce the functions 



A?(t) = ±J**?(s)d8, i = A, BA, 



d2 



which together form a 2D vector A N (t) . We will also denote by Nq^(t) G 
the pair of queues at the moment t. By NQf(t), i = A,BA we denote the 
number of clients which have arrived to the server A during the time interval 
[0, t] . We can as well assume that to every client the service time is assigned 
at the moment of its arrival. The sum of these required service times for 
clients arrived during the time interval [0, t] will be denoted by NW t N (t); this 
function has a stair-like graph. By Nwf(t) we denote the remaining required 
times for clients queuing or being served at the moment t; the graph of these 
functions are saw-like. 

We also consider the fluid model at this node. So (t) will denote the 
evolution of the initial measure q^ (0) on under fluid dynamics governed 
by the inflow with rate A^ (t) . By Qi{t) = A«(t) we denote the amounts of 
fluid arriving to our node during the time interval [0, t\. Again, the fluid BA 
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has priority over the fluid A, that is, it goes out first at rate 71 whenever 
present. 

Our goal is to prove convergence to the fluid limit. Let us assume that 
the scaled inflows and initial states of the node converge to those of the fluid 
model as N — * 00, that is 

lim sup || A* (t) - A(t)|| = (41) 

Ar - > °° te[o,T] 

and 

lim pkrov (q N {0),q(0)) = 0. (42) 

In what follows, we will need the following three bounds. The first one is 
the statement that the process W N (t) is very close to the function r y~ 1 A N (t) 
- namely, 

lim E [ sup \\W N (t)- 1 ' 1 A N (t)\\ j = 0. (43) 



te[o,T] 



A simpler claim concerns the sum Nw N (0) of the service times of all the 
users present in the queue Nq N (0) at the initial moment t = 0. Namely, for 
the conditional distribution of w N (0) under the condition q N (0) = q we have 



E 



u^(0)|/ r (0) = g )- 7 - 1 g| <m)h\\ (44) 



for some ip(N) — > as N — > 00. Moreover, 

sup pkrov {w N (t), 1 - l q N (t)) < ip(N) sup pkrov (q N (t),0) . (45) 
*e[o,r] te[o,T] 



Of course, for every fixed t the convergence in (J43J) follows from the Central 
Limit Theorem. The problem is that we need the convergence at all moments t. 
We will obtain (1431) by constructing the finite-point event, which contain the one 
we are interested. Indeed, consider the event E (t , c) , c > 0, which consists of 
all trajectories such that 

(to) > Ta 1 (A£ (t ) + c) . (46) 

Note that the function A^ (t) has its derivative < C (see (I40p ). so if the event 
E (to, c) happens, then all the events E (t + r, c — Ct) happen as well, since the 
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function W% (t) is non-decreasing on every trajectory. Therefore we can replace 
the infinite union by the finite one: 



U E(t ,c)C 

0<t <T 



2CT/c 



{)E(k 

k=0 



c c 
2C , 2 



and use the fact that for any fixed t < T and c the probability of the event 
E (t, |J is exponentially small in N. The remaining cases (corresponding to the 
second coordinate, AB, and to the lower estimate in (146 P ) are immediate. 

The proof of (1451) proceeds in a similar way. First of all, there is a natural 
coupling between the processes w N (t) and q N (t), corresponding to the fact that 
we can assume that the service time of every client is known at its arrival moment. 
Consider the event E' (to, c) , consisting of the trajectories where 

<(*>)> 7a 1 (<#(*>) + 

We would like to use the argument similar to the above, saying that if the bad 
event E' (to, c) happens at to, then on a whole segment around to something 
unlikely has to happen as well. Partially it can be done, since for every tra- 
jectory we have w^(to + t) > w^(t ) — r, and so we have for every trajec- 
tory in E' (t , c) that w^(t + r) > 7^ ((£4 (t ) + c — 7A r ) • If we can claim 
that (/a (to) + c — 7at > (/a (to + t) + c/2 for all r small enough, then we 
would be done. However, the outcome that q^(to + r) > q^(to) + c /2 — 
7at is not excluded, even if r is very small. Yet, the probability of the in- 
crease of the queue by c/2 — 7^4 t during the time r is exponentially small in 
r as t — > 0. Therefore, the event Uo<t <T-E" (to, c) is contained in the union 



U 



2C'T/c 
k=0 



R' (k— -) 

^ V^ iC'i 2) 



U 



u 



2C'T/c j 
fc=0 



[k-^jj) f ) for some suitably chosen C, where 



E" [k-^y-, |) is the event that on the segment [k-^h, (k + 1) ^r] the increment 

?A ((^ + 1) ~~ (k^t 7 ) > c / 4 ' and we are done - 
From OH]) and (jig]) , we get 



lim E sup llW^t) -7" 1 A(t)|| = 0. 



JV— »oo 



te[o,T] 



Note that 

taov K(0), 7" 1 ?(0)) < E (E I (w N (0) I g*(0) = q 

+ PKROV (7" 1 g 7V (0),7" 1 g(0)) 



(47) 



7 1 q 
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(Here we treat w N (0) in the l.h.s. as the probability distribution.) So we get 
from (HZD and (B3D that 



Pkrov K(0), 7 _1 g(0)) < <p(N)p KROV (g(0), 5„) , (48) 

where y(iV) — > as AT — ► oo, and GM 2 is the origin. 

Next, we need the following estimate (see for instance [LTj ) : 

Lemma 13 Let A(t) and A'(t) be two inflows to the fluid priority node with 
initial (non-random) fluid levels q(0) and q'(0). Then 

sup \\ q ( t )- q '(t)\\<L(\\ q (0)-q'(0)\\+ sup ||A(t)-A'(t)||), 
tG[0,T] te[o,T] 

where L = L(7i,7 2 ). 

Proof. Let us consider a fluid single-class node with variable capacity. 
Namely, let ?i(0) be the (scalar) initial fluid level, let Ai(t) be the inflow, and 
introduce S\(t) to be the server capacity, which is the amount of work the 
server can do during the time interval [0, t\. (For example, in our situation 
S\{t) = lit, but we will consider more general case, with Si(t) not necessarily 
linear.) Introduce the virtual level 

V{t) = q x {0)+A 1 {t)-3 1 {t) (49) 

and the unused service capacity 

U(t) = max{0, - inf Vis)}. (50) 
se[o,i] 

Then 

qi(t) = V(t) + U(t), t>0. (51) 

Let us introduce the sup norm on the space of functions. Then the functionals 
{ gi (0), AxG), ^ V(-) and {V(-), U(-)} ^ Ql (-), given by fl3B and (JEEJ, 

have finite norms, since they are linear. The non-linear functional V(-) — > 
U(-), given by fl50|) . has finite norm as well. Indeed, the functional V(-) — > 
inf[ .] V(-) has norm < 1, since for any pair x(-),y(-) of scalar functions 



inf x(s) — inf -u(s' 
«e[o,t] se[o,t] 



< sup |x(s) — y(s) 
«e[o,t] 
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and so 



sup 

te[o,T] 



inf x(s) — inf y(s) 
se[o,t] se[o,t] 



< sup |x(t) - y(t)|; 
te[o,T] 



the same holds for the functional {x (t)} —>■ max{0, x(t)}, since 

sup |max{0, x(t)} — max{0, y(t)}\ < sup \x(t)—y(t)\. 

te[o,r] *e[o,T] 

Therefore the composed functional, taking the triplet {?i(0), Ai(-), Si(-)} to 
the pair {qi(-), U(-)}, defined by f|49|) - fl5T]) . has finite norm. 

That proves the desired statement for the first component of q. Now, to 
finish the proof, we note that the capacity of the server for the users of the 
second class is given by ^(t) = j2U(t), where U(t) is the unused server 
capacity for the high-priority class (with S±(t) = jit, t > 0). Then we repeat 
the argument above. ■ 

Next we formulate as a separate statement the obvious remark that the 
evolution of the current remaining service time variable, w (t) , coincides with 
the evolution of the level of some evidently constructed fluid system. 

Lemma 14 Let users Uj of two possible types i = A, BA, with service times 
h\ arrive at the initially empty server at times t\, j — 1, 2, . . . , and let Wi(t) 
be the evolution of the remaining service times. Consider also the fluid model 
with two classes of fluids, which starts in the empty state and is governed by 
the fluid inflows 

Ai(*) = 7< Yl h i- 

(I.e., our fluids have "viscosities" 7f : and^ 1 -) Then at every momentt > 
the current levels of fluids at the server equal to (i), i = 1,2. 

Another auxiliary result is needed: 

Lemma 15 Let q G 1R 2 be (random) queue to our server, and w be the 
corresponding (random) amount of total work (=service time needed). The 
service times of the users are independent, and within the i-th class identically 
distributed with mean 7" 1 . Then, for any v E R 2 , 

Pkrov {l^q, 5v) < Pkrov (w, 5 v ) . (52) 
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Proof. Since the norm II • II is convex, we have for the conditional random 



variable w q = q that ||7 1 q — v\\ < E|| (^w 

gives ([52]) . ■ 

Now, we derive the following result 

Proposition 16 Under the assumptions made above, 



q = q) — v\\. Averaging over q 



lim sup Pkrov (q N {t) , = 0- (53) 



Proof. First, 



sup Pkrov {q N (t),q(t)) <max{7i,7 2 } sup Pkrov (l 1 q N (t),l 1 q(t)) . 
*€[0,T] te[o,T] 

Then we write a chain of inequalities. First of all we have 

sup pkrov ( r Y~ 1 Q N (t)n~ l q(t)) 
te[o,T] 

< sup pkrov (w N (t),^ 1 q N (t)) + sup p KROV (w N (t),^ 1 q(t)) , 
te[o,T] te[o,T] 

and the first summand can be bounded by ( ]4"5i) . To estimate the distance 
Pkrov {w N (t), 7 we have to exhibit some joint distribution of 7~ 1 g(t) 
and w N (t). We take the following one: first, we choose the coupling between 
7 _1 g(0) and w N (0), using (jUj) and (H2D , getting 

Pkrov (w n (0), 7 _1 ?(0)) < 2^(N) Pk rov (q(0), 0) . 

Given the joint realization of the initial values (7~ 1 q , (0), w N (0)) , the evolu- 
tion of the coordinate q(t) is deterministic, defined by the flows A^(i) of the 
arriving fluids. The evolution w N (t) is stochastic, governed by the Poisson 
process with net rates (i). Therefore for every t we have to exhibit the 
coupling between the distribution of the vector w N (t) and deterministic value 
q(t). The resulting KROV distance is precisely what the Lemmas [TBI and [T 
allow us to control: 



PKROV [w"^-^® w N (0),q(0)^ 
= E(||«; JV (t)-7- 1 g(*)|| \w N (0),q(0)) 



< L ( 11^(0) - 7" 1 g(0)|| + E ( sup \\W N (s) - 7 _1 A(s) 

se[o,t] 
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It remains to apply bounds (H7|) and (14"5|) and deduce (15^1) . ■ 

2. The proof of our statement for the closed system does not require any 
extra arguments, since the closed system is a special case of the open system, 
where all the flows satisfy the relations defining the closed system. 



6 Main result 

In this section we finally formulate and prove our main theorem, which claims 
that the NLMP started from some special initial state behaves similarly to 
the fluid model in its periodic regime, at all times t G (0, oo). Since we 
know already that the NLMP is in turn a limit of networks of size M, as 
M — > oo, our theorem implies that the large size (M ^> 1) Markov process 
behaves similarly to the fluid model for a very long time, which time 
diverges as M — > oo, provided the number N of clients per node exceeds some 
value N . In particular, there are initial states for the networks Vj^-, which 
lead to a long time oscillations, before the network reaches its stationary 
state. 

Theorem 17 Let e > 0. Then there exist the values N , e' > 0, a > and 
E < oo such that for all N > N the states v N (t) of the NLMP process V£, 
started at the initial state v N (0) with the properties: 

Pkrov {y N (0) A(o)) < ef, 

with x (0) G C, 

(exp{aL(x)})^ (0) < E, 

satisfies for allt > 

Pkrov (t) , S Xv ) < e, (54) 

where x u G C is some moving point, depending on the process v = lv N (t) , t > 
In particular, the process u N (t) has no limit as t —>■ oo. 

In words, we are proving that if we start the NLMP with high load iV 
per server, from the state close to some atomic measure 5 Z with z belonging 
to the cycle, then it never goes to a limit. 

For that, we need a general Lemma, which is formulated in the Euler 
scaling. First, we recall the definitions. Let A (t) = {Xi (t) ,i — 1, k} be 
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the rates of Poisson inflows of the customers of k types, and 7$, i — 1, k 
be their rates of service. The discipline of service will be irrelevant here; we 
need only that the server is not idle if the queue is not empty. We call the 
flow A (t) to be underloaded, with parameters (T, 5) , if for any t 



i=l ll Jt 



t+T 

K (t) dt<(l- 5) T. 



Below we are talking about the flow with load N. That means that we con- 
sider the situation when the input rates are given by NX (t) = {NXi (t) , i — 1, . 
while the service rates equals to N^. 

Lemma 18 Consider the Non-Homogeneous Markov Process (J,(t), started 
from the initial state /J, (0) , and suppose that its generating rate function 
A (t) is underloaded, with parameters (T, 5) . Then there exist values a > 
and A < oo, depending only on the pair (T, 5) , such that if the exponential 
moment (exp {aL (x)}) of the initial state is finite, then for times t > 
t (/j, (0)) and for any load N 

(exp{aL(x)}) A4(t) < A. 

Moreover, there exists the time T = T (T, 5) , such that for any initial state 
/i(0) 7 satisfying the estimate 

(exp{aL(x)}) m < 3 A, 

we have, for any load N, that 

(exp{aL(x)» MT) < 2A. 

Proof. We consider Poisson inflow with a general distribution rj of the service 
time, having finite exponential moment. In particular, exponential service 
time fits. 

The users arrive to the node according to the Poisson processes with 
rates Aj(t). Their service times are i.i.d. with distribution functions f]i(h). 
Let Ei = E (r/i) . We study the dynamics of the remaining service time, hence, 
the service discipline is of no importance. The regime we are interested in is 
the underloaded regime; that means that for some 5 > 0, all T large enough 
and all t > 

f t + (j2Hs)E t )ds<T(l-5). 
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Let dfi t {u) be the current distribution of the remaining service time. We 
want to study the exponential moment 



Q a {t) = / q a {u)dnt{u), (55) 
Jo 

where q a (u) = e au . We will show that the moment Q a (t) satisfies the equa- 
tion 

Q a {t) < e Cl -^Q a (0) + C 2 , t>0. 

The statistics of the observable L will then be easy to derive. 

Note that the underload condition ensures the absolute continuity of 
Q a (t). We will need the quantities 



POO 

$f= / (e ah - 1) drji(h). 
Jo 



We assume that < +00 for all a < a with a > 0. 

Before studying the moments (1551) , we will consider the situation of the 
"broken" server, when the clients (of one type) only come, but are not served. 
The queue then only grows in time, as is the workload u. The corresponding 
exponential moment will be denoted by Q a if). We have: 

Q^(t) = X(t)^Q^(t). (56) 

Indeed, the event of arrival of a user with service time h at the queue with 
the current workload u shifts the workload to the value u + h, so the value of 
q a changes from e au to e a ^ + ^ = e au + e au ( e ah — 1). In order to find Qa (t), 
we have to multiply the increment e au (e ah — 1) by the rate X(t) of the arrival 
event and to integrate it with respect to dfi t {u) x dr](h), since u and h are 
independent. In this way we arrive to fl56|) . 

Next, let us study the case of "broken pipe", when the inflow is zero, so 
the server works only on the initial supply of clients. The evolution of the 
distribution fi t of the workload is given by the following simple relation: 



Vt+s [a, b] 



fit [a + s, b + s] if < a < b, 
yU t (— 00, b + s] if a < < b. 



In words, the atom at u = grows with time. We denote the corresponding 
exponential moment by Qa (t). The straightforward computation shows that 

Q^(t)=a Po (t)-aQ^(t), (57) 
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where Po(t) is the current probability of the queue to be empty. Note that 
Qa\t) is absolutely continuous and flBTI) holds for almost all t. 

In the general case of several inflows we put the two relations together to 

get 



Qa(t) 



a 



Q a {t) + ap Q (t). 



(58) 



Let us now rewrite $?. We have 



POO 

= / {e ah - 1) d Vi (h) = 
Jo 

poo poo 

/ ah dr)i(h) + / [e ah - 1 - ah] drn(h) = aEi + aF { \ 
Jo Jo 



a) 



(59) 



where E± is the mean service time and F^a) is continuous function of a G 
[0, a], which satisfies Fj(a) = O (a) as a — > 0. From f[55|) and ([59]) we get for 
a < 1 the bound 



&*(*) < a 



a) 



Q a {t) + 1. (60) 



The Euler scaling with parameter N changes A(i) to A iV (t) = NX(t) and 
^(/i) to r/ N \h) = Nrj(Nh). Hence, Xf(t)E^ does not depend on N. Let us 
show that NFf (a) is small for all N, once a is small. Indeed, 

POO 

NaF N (a) = N 2 [e ah -I- ah] dr](Nh) 
Jo 



N 2 



e §Nh _ 1 



a 



Nh 



drj(Nh) 



Hence, we get a uniform bound for a small enough and all iV > 1 simul- 
taneously (and, by the limit, for the fluid model U N = oo" as well): 



QZ(t)<a 



J2 X i(t) {Ei + K (a))-l 



(61) 



where x (a) ~ a. 
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The solution to the linear equation 



x(t) 



a{t)x(t) + b (t) 



is given by the formula 



x(t) 



g(0,t)x(0) + [ g(s,t)b(s)ds, 



where g(s,t) = e^ a{r)dT . We apply it to (JBTJ), with x (t) = Q a (t), a(t) = 
a K(t) (Ei + x (a)) — 1] and b(t) = 1. By the underload assumption, 



In our case the distribution r\i is exponential with the parameter 7^. Let 
us show finally that the exponential bound on the workload implies an expo- 
nential bound on the number of customers (may be, with another exponent). 

Indeed, under the condition that we are in the state with rii and n<i 
customers of two classes, the conditional distribution of the workload u is a 
measure [i ni n 2 011 , with mean value u = ^ + By convexity of the 

exponent, 




for some C\, j3 > and for all s < t, once a is small. Then, J Q g(s, t)ds < C2 
for all t > 0. 
Hence, 

Q a (t) < e Cl -^Q Q (0) + C 2 , t>0. 




au 



which provides us with the upper bound 




Taking expectations with respect to n\ (t) , 112 (t) we get 



E ( eTl+T2 



(m(t)+n 2 (t)) 



) 



< 



/ 



e au dfi t (u) = Q a (t), 



which is the desired estimate. 
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Proof of the Main Theorem. The proof proceeds by "induction" in 
time. We suppose inductively that at a certain (Euler) moment T the NLMP 
with the load N is in the state fi (T) , having two properties: 

(exp{aL(x)}) KT) <3A, (62) 

Pkrov (p (T) , S x ) < e (63) 

for some x G C. We will show that there exists the time T', at which the 
same two conditions hold for the measure p. (T + T') - except for different 
point x on the cycle C. 

To see this we first consider the Non-Linear Dynamical System (NLDS) 
Aqo, with initial state fj, (T) . In other words, Aoo = A^ (Y (■)) , for some 
Y (•) G y(fj,(T)). We can use the Proposition [TUl which tells us that for 
any T' large enough pkrov (A^/x (T) ,5 Z (T')) < £/3, Choosing one such T" 
(uniformly in /j (T) , satisfying (I6"2l — (1631 !) we can claim that for the NLMP 
evolution we have Pkrov {p (T + T') , 5 x (t')) < 2e/3, provided only that 
is large enough; indeed, we know from Theorem [121 that the NLMP converges 
in the KROV metric to NLDS on any finite time interval, as iV — > oo, which 
convergence is uniform over the set of initial measures satisfying (|62|) . Note 
that we thus have reproduced the condition (|63l) . 

Proposition [TT] tells us that for all t < T' pkrov 

(A^/i (T) , 6 x (t)) < 

e(T',e) . Due to the same convergence statement, for the NLMP evolution 
we have Pkrov (p (T + 1) , 5 z (t)) < 2e(T', e) for all t < T'. In words, the 
measure /jl(T + t) goes very close to the cycle trajectory. Since we want 
to use the Lemma [181 we can as we ll assume that T' > T, where T is 
the time introduced in this Lemma. Now all its conditions are satisfied, so 
Lemma [TBI tells us that (exp{aL(x)}) , T+T n < 2A, thus the condition f[6"2l 
is reproduced as well. ■ 

7 Conclusions 

Our main result indicates that there is an important similarity between large 
queuing networks and large systems of statistical mechanics. Namely, we 
have shown that the load per server plays for some networks the same role as 
the inverse temperature in statistical mechanics. At high load the network 
can lose the property of uniqueness of the stationary state and start to behave 
in the oscillatory manner. This phenomenon looks similar to the fact that 



50 



some 3D systems with continuous symmetry are not ergodic under Glauber 
dynamics, when the temperature is low enough. 

It is very interesting to understand how general this phenomenon is; our 
expectations are that such non-ergodic behavior is a characteristic feature of 
the high load regime. 

In the forthcoming publications we will show that the behavior in the low 
load regime is always ergodic, which corresponds to the high temperature 
uniqueness of statistical mechanics. 
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